###MC experiment of Lin 2024

rm(list = ls())
ptm <- proc.time()

setwd('.../Monte Carlo')

# library(foreach)
# library(doParallel)
# numCores <- detectCores()
# registerDoParallel(numCores) 
set.seed(12345);
n=400;
loop=1000;
tau=0.5


sms=matrix(rep(0,5*loop),loop,5);

delta=-1;
##True parameter     
beta=c(1,1,1,1);
beta=as.vector(beta);
alpha=1;
gamma=1;

library(Matrix)




# sms= foreach(l=1:loop, .combine=rbind) %dopar% {
for(l in 1:loop){   cat(l,"Iteration \n")
  ##Generating Large network with n
  FF=matrix(rep(0,10*10*n),10*n,10);
  for(g in 1:n){
    NF = sample(0:10,10,TRUE);##Generating Friendship with maximum 10 friend.
    for(i in 1:10){
      location=sample(1:10,NF[i],TRUE);
      FF[10*(g-1)+i,location]=1;
      FF[10*(g-1)+i,i]=0;
    } 
  }
  
  Friend=FF[1:10,]
  for(g in 1:(n-1)){
    Friend=bdiag(Friend,FF[(10*g+1):(10*g+10),])
  }
  
  Friend=as.matrix(Friend)
  NF = rowSums(Friend);
  invNF =1/NF;#1 over number of friend
  invNF[which(!is.finite(invNF))]=0;
  
  x1 = rnorm(10*n);
  x2 = 2*rbinom(10*n,1,0.5)-1;
  x3 = rt(10*n,5);
  xx=cbind(1,x1,x2,x3);
  x=cbind(x1,x2,x3);
  
  u=rlogis(10*n,scale=1);
  
  ##Find the true BNE with fixed point algorithm
  t=10
  P0=rep(0,10*n);
  P1=rep(0,10*n);
  tot=100;
  reps=0;
  while (t > 1E-2&reps<tot){
    reps=reps+1;
    P1 = as.integer(xx%*%beta+alpha*as.matrix(Friend)%*%as.vector(P0)*as.vector(invNF)>0);
    t=mean(abs(P1-P0))
    P0=P1;
  }
  #Latent true decision
  y = as.integer(xx%*%beta+alpha*as.matrix(Friend)%*%as.vector(P1)*as.vector(invNF)-0.5*u>0);
  m=mean(y)
  
  bne=rep(0,10*n);
  xnet0=as.matrix(Friend)%*%as.matrix(bne)*invNF; ##Construct network regressor
  
  
  source("sa.nps.R")
  err=10;
  err.tol=1E-2;
  bne0=rep(0,10*n);
  npsms.reps = 0;
  tot.reps = 100;
  quiet=FALSE;
  nc=0;
  while(err>err.tol&npsms.reps<tot.reps){
    npsms.reps = npsms.reps + 1; 
    npsms=sa.nps(y,x,xnet0,tau,sigma=0.5*(10*n)^(-1/5),lbc=-10,ubc=10);

    betahat0=npsms$coef;
    xx1=cbind(1,xnet0,x);
    bne1 = as.integer(as.matrix(xx1)%*%as.vector(betahat0)>0);
    err = mean(abs(bne1-bne0))
    xnet0=as.matrix(Friend)%*%as.matrix(bne1)*invNF;
    bne0=bne1;
  }
  sms[l,]=as.vector(t(npsms$coef))
}



###AB
ABa=mean(sms[,2])-alpha
ABb=colMeans(sms[,c(1,3:4)])-beta[1:3]

###MSE
MSEa=sum((sms[,2]-alpha)^2)/loop

MSEb=colMeans((sms[,c(1,3:4)]-t(matrix(rep(beta[1:3],loop),3,loop)))^2)


round(c(ABb,ABa),digits=3)
round(c(MSEb,MSEa),digits=3)


proc.time() - ptm
